<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>OptForce GAMS</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S4 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S5 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }
.S6 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S16 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>OptForce GAMS</span></h1><h2  class = 'S1'><span>Author: Sebastián N. Mendoza,  Center for Mathematical Modeling, University of Chile. snmendoz@uc.cl</span></h2><h2  class = 'S1'><span style=' font-weight: bold;'>Reviewers(s): Chiam Yu Ng (Costas D. Maranas group), </span><span>Lin Wang </span><span style=' font-weight: bold;'>(Costas D. Maranas group)</span></h2><h2  class = 'S1'><span style=' font-weight: bold;'>INTRODUCTION:</span></h2><div  class = 'S2'><span>In this tutorial we will run optForce. For a detailed description of the procedure, please see [1]. Briefly, the problem is to find a set of interventions of size "K" such that when these interventions are applied to a wild-type strain, the mutant created will produce a particular target of interest in a higher rate than the wild-type strain. The interventions could be knockouts (lead to zero the flux for a particular reaction), upregulations (increase the flux for a particular reaction) and downregulations (decrease the flux for a particular reaction). </span></div><div  class = 'S2'><span>For example, imagine that we would like to increase the production of succinate in Escherichia coli. Which are the interventions needed to increase the production of succinate? We will approach this problem in this tutorial and we will see how each of the steps of OptForce are solved. </span></div><h2  class = 'S3'><span>MATERIALS</span></h2><h2  class = 'S1'><span>EQUIPMENT</span></h2><ol  class = 'S4'><li  class = 'S5'><span>MATLAB</span></li><li  class = 'S5'><span>GAMS</span></li><li  class = 'S5'><span>A solver for Mixed Integer Linear Programming (MILP) problems. For example, Gurobi.</span></li></ol><h2  class = 'S1'><span style=' font-weight: bold;'>EQUIPMENT SETUP</span></h2><div  class = 'S2'><span style=' font-weight: bold;'>GAMS</span><span>: Install the latest version of GAMS available. This is important as some GAMS functions are not supported in versions older than  24.7.1. During GAMS installation make sure to select the option "Add GAMS directory to PATH environment variable". In MATLAB, use functions addpath() and savepath() for adding and saving the directory where GAMS was installed. For a more detailed description of GAMS installation and setup, you can watch the official tutorial of GDXMRW at https://www.youtube.com/watch?v=HPn_q8nlktE</span></div><h2  class = 'S3'><span>PROCEDURE</span></h2><div  class = 'S2'><span>The proceduce consists on the following steps</span></div><div  class = 'S2'><span>1) Maximize specific growth rate and product formation.</span></div><div  class = 'S2'><span>2) Define constraints for both wild-type and mutant strain: </span></div><div  class = 'S2'><span>3) Perform flux variability analysis for both wild-type and mutant strain.</span></div><div  class = 'S2'><span>4) Find must sets,  i.e, reactions that MUST increase or decrease their flux in order to achieve the phenotype in the mutant strain. </span></div><h2  class = 'S1'><span>Figure 1.</span></h2><div  class = 'S6'><img class = "imageNode" src = "" width = "404.1" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>5) Find the interventions needed that will ensure a increased production of the target of interest</span></div><div  class = 'S2'><span>Now, we will approach each step in detail.</span></div><h2  class = 'S3'><span>STEP 1: Maximize specific growth rate and product formation</span></h2><div  class = 'S2'><span>First, we load the model. This model comprises only 90 reactions, which describe the central metabolism of E. coli [2].</span></div><div  class = 'S2'><span>Then, we change the objective function to maximize biomass ("R75"). We also change the lower bounds, so E. coli will be able to consume glucose, oxygen, sulfate, ammomium, citrate and glycerol.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">global </span><span >TUTORIAL_INIT_CB;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~isempty(TUTORIAL_INIT_CB) &amp;&amp; TUTORIAL_INIT_CB==1</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'gurobi'</span><span >,</span><span style="color: rgb(170, 4, 249);">'all'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >pathTutorial = which(</span><span style="color: rgb(170, 4, 249);">'tutorial_OptForceGAMS.mlx'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >pathstr = fileparts(pathTutorial);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >cd(pathstr)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'AntCore.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >modelDirectory = getDistributedModelFolder(modelFileName); </span><span style="color: rgb(2, 128, 9);">%Look up the folder for the distributed Models.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >modelFileName= [modelDirectory filesep modelFileName]; </span><span style="color: rgb(2, 128, 9);">% Get the full path. Necessary to be sure, that the right model is loaded</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = readCbModel(modelFileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model.c(strcmp(model.rxns,</span><span style="color: rgb(170, 4, 249);">'R75'</span><span >)) = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_gluc'</span><span >, -100, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_o2'</span><span >, -100, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_so4'</span><span >, -100, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_nh3'</span><span >, -100, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_cit'</span><span >, -100, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_glyc'</span><span >, -100, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span></span></div></div></div><div  class = 'S10'><span>Then, we calculate the maximum specific growth rate and the maximum production rate for succinate</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >growthRate = optimizeCbModel(model); </span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'The maximum growth rate is %1.2f'</span><span >, growthRate.f);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="985FE36B" data-testid="output_0" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">The maximum growth rate is 14.36</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = changeObjective(model, </span><span style="color: rgb(170, 4, 249);">'EX_suc'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >maxSucc = optimizeCbModel(model);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'The maximum production rate of succinate is %1.2f'</span><span >, maxSucc.f);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="FBFA982C" data-testid="output_1" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">The maximum production rate of succinate is 155.56</div></div></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>TIP: </span><span>The biomass reaction is usually set to 1%-10% of maximum theoretical biomass yield when running the following steps, to prevent solutions with not biomass formation</span></div><ol  class = 'S4'><li  class = 'S5'><span>maximizing product formation</span></li><li  class = 'S5'><span>finding MUST sets of second order</span></li><li  class = 'S5'><span>finding FORCE sets</span></li></ol><h2  class = 'S3'><span>STEP 2: Define constraints for both wild-type and mutant strain</span></h2><div  class = 'S2'><span style=' font-weight: bold;'>TIMING</span><span>: This step should take a few days or weeks, depending on the information available for your species. </span></div><div  class = 'S2'><span style=' font-weight: bold;'>CRITICAL STEP</span><span>: This is a manual task, so you should search for information in articles or even perform your own experiments. You can also make assumptions for describing the phenotypes of both strains which will make this task a little faster but make sure to have two strains different enough, because you should be able to find differences in reactions ranges. </span></div><div  class = 'S2'><span>First, we load the model. This model comprises only 90 reactions, which describe the central metabolism of E. coli [2].</span></div><div  class = 'S2'><span>Then, we change the objective function to maximize biomass ("R75"). We also change the lower bounds, so E. coli will be able to consume glucose, oxygen, sulfate, ammomium, citrate and glycerol. </span></div><div  class = 'S2'><span>We define constraints for each strain</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >Constr_WT = struct(</span><span style="color: rgb(170, 4, 249);">'rxnList'</span><span >, {{</span><span style="color: rgb(170, 4, 249);">'R75'</span><span >}}, </span><span style="color: rgb(170, 4, 249);">'rxnValues'</span><span >, 14, </span><span style="color: rgb(170, 4, 249);">'rxnBoundType'</span><span >, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >Constr_MT = struct(</span><span style="color: rgb(170, 4, 249);">'rxnList'</span><span >, {{</span><span style="color: rgb(170, 4, 249);">'R75'</span><span >,</span><span style="color: rgb(170, 4, 249);">'EX_suc'</span><span >}}, </span><span style="color: rgb(170, 4, 249);">'rxnValues'</span><span >, [0, 155.55], </span><span style="color: rgb(170, 4, 249);">'rxnBoundType'</span><span >, </span><span style="color: rgb(170, 4, 249);">'bb'</span><span >);</span></span></div></div></div><h2  class = 'S3'><span>Step 3: Flux Variability Analysis</span></h2><div  class = 'S2'><span style=' font-weight: bold;'>TIMING</span><span>: This task should take from a few seconds to a few hours depending on the size of your reconstruction</span></div><div  class = 'S2'><span>We  run the FVA analysis for both strains</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[minFluxes_WT, maxFluxes_WT, minFluxes_MT, maxFluxes_MT,~,~] = FVAOptForce(model, Constr_WT, Constr_MT);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: normal"><span >disp([minFluxes_WT, maxFluxes_WT, minFluxes_MT, maxFluxes_MT]);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="8A5D6CFD" data-testid="output_2" data-width="428" data-height="1263" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  -90.1251   97.1300   44.4313  100.0000
         0   86.0700   44.4375  100.0000
         0   86.0700   44.4375  100.0000
  -56.1567   86.0700  -44.4500   11.1143
   21.3033  163.5300   55.5500  111.1143
   -3.0777  154.8640   55.5500  111.1143
         0  151.5086         0   55.5625
         0  187.2551         0   55.5687
         0  169.5163         0    0.0187
  -10.0660  102.9449         0    0.0125
   10.0660   66.5714         0    0.0063
  -10.0660  102.9449         0    0.0125
  -48.9454    7.5600   -0.0063         0
  -53.9994    2.5060   -0.0063         0
  -53.9994    2.5060   -0.0063         0
   -2.5060   53.9994         0    0.0063
         0   86.0700         0   55.5625
         0   86.0700         0   55.5625
    9.7020  114.6466   55.5500   55.5625
         0   56.5564   55.5500   55.5571
   16.0264  145.2048  155.5500  155.5563
   16.0264  145.2048  155.5500  155.5563
    0.9344  130.1128  155.5500  155.5562
   -5.6736  123.5048  155.5500  155.5563
         0  118.0576         0    0.0062
    5.1940  123.2516         0    0.0062
  -98.1150  123.2516  -55.5625    0.0062
         0  151.5086         0   55.5625
         0  151.5086         0   55.5625
         0  254.5400   55.5500  777.7875
         0  253.2493         0  722.2375
   -7.1960   94.6056         0    0.0125
         0   84.8467   88.8750   88.9000
         0   84.8467   88.8750   88.9000
         0  175.1064  188.8500  188.9000
         0  175.1064  188.8500  188.9000
   91.4130  107.1280         0         0
    9.4500    9.4500         0         0
    2.9400    2.9400         0         0
    3.9340    3.9340         0         0
   25.4520   56.8820         0         0
    3.2060    3.2060         0         0
    6.8320    6.8320         0         0
         0   15.7150         0         0
   -6.8880    8.8270         0         0
    0.6790   16.3940         0         0
         0   31.4300         0         0
    3.2620    3.2620         0         0
    4.5640    4.5640         0         0
    4.5640    4.5640         0         0
    7.2380   38.6680         0         0
    2.0440    2.0440         0         0
    5.6280    5.6280         0         0
    5.9920    5.9920         0         0
    3.8640    3.8640         0         0
    2.4640    2.4640         0         0
    1.8340    1.8340         0         0
    0.7560    0.7560         0         0
    1.2600    1.2600         0         0
    2.0440    2.0440         0         0
    1.2600    1.2600         0         0
   79.7324  200.0000  199.9500  200.0000
         0  118.0576         0    0.0062
  -39.5563  353.9124  -22.2500   33.3500
         0  253.2493         0  722.2375
   40.6268  100.0000   99.9875  100.0000
   15.0890  100.0000   99.9929  100.0000
 -100.0000   84.8467 -100.0000  -99.9500
         0  175.1064  188.8500  188.9000
         0  101.8016         0    0.0125
  134.9718  407.3274  311.1000  311.1187
   62.1267  100.0000   99.9750  100.0000
   97.4820   97.4820         0         0
    3.2620    3.2620         0         0
   14.0000   14.0000         0         0
         0  175.1064  188.8500  188.9000
  134.9718  407.3274  311.1000  311.1187
         0  101.8016         0    0.0125
         0  253.2493         0  722.2375
 -100.0000  -40.6268 -100.0000  -99.9875
 -100.0000  -15.0890 -100.0000  -99.9929
 -100.0000   84.8467 -100.0000  -99.9500
  -97.4820  -97.4820         0         0
 -100.0000  -62.1267 -100.0000  -99.9750
   -3.2620   -3.2620         0         0
         0  105.4230  155.5500  155.5500
         0  105.4230  155.5500  155.5500
   11.6200   11.6200         0         0
    5.0540    5.0540         0         0
    5.9920    5.9920         0         0</div></div></div></div></div><div  class = 'S10'><span>Now, the run the second step of OptForce.</span></div><h2  class = 'S3'><span>Step 4: Find Must Sets</span></h2><div  class = 'S2'><span style=' font-weight: bold;'>TIMING: </span><span>This task should take from a few seconds to a few hours depending on the size of your reconstruction</span></div><div  class = 'S2'><span>First, we define an ID for this run. Each time you run the functions associated to the optForce procedure, some folders can be generated to store inputs used in that run. Outputs are stored as well. These folder will be located inside the folder defined by your run ID. Thus, if your runID is ''TestOptForce", the structure of the folders will be the following:</span></div><div  class = 'S2'><span style=' font-family: monospace;'>├── CurrentFolder</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   ├── TestOptForce</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── Inputs</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   └── Outputs</span></div><div  class = 'S2'><span>To avoid the generation of inputs and outputs folders, set keepInputs = 0, printExcel = 0, printText = 0 and keepGamsOutputs = 0</span></div><div  class = 'S2'><span>Also, a report of the run is generated each time you run the functions associated to the optForce procedure. So, the idea is to give a different runID each time you run the functions, so you will be able to see the report (inputs used, outputs generated, errors in the run) for each run.</span></div><div  class = 'S2'><span>We define then our runID</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: normal"><span >runID = </span><span style="color: rgb(170, 4, 249);">'TestOptForce'</span><span >;</span></span></div></div></div><div  class = 'S10'><span>Fow now, only functions to find first and second order must sets are supported. As depicted in </span><span style=' font-weight: bold;'>Figure 1</span><span>, the first order must sets are MUSTU and MUSTL; and second order must sets are MUSTUU, MUSTLL and MUSTUL</span></div><div  class = 'S2'><span style=' font-weight: bold;'>A) Finding first order must sets</span></div><div  class = 'S2'><span>We define constraints</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: normal"><span >constrOpt = struct(</span><span style="color: rgb(170, 4, 249);">'rxnList'</span><span >, {{</span><span style="color: rgb(170, 4, 249);">'EX_gluc'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R75'</span><span >,</span><span style="color: rgb(170, 4, 249);">'EX_suc'</span><span >}}, </span><span style="color: rgb(170, 4, 249);">'values'</span><span >, [-100, 0, 155.5]');</span></span></div></div></div><div  class = 'S10'><span>We then run the functions findMustLWithGAMS.m and findMustUWithGAMS.m that will find mustU and mustL sets, respectively.</span></div><div  class = 'S2'><span>Important: To run these function you will need a solver able to solve Mixed Integer Linear Programming (MILP or MIP) problems. Some popular options are: cplex and gurobi. You can see which gams solvers are available in your systems to solve MIP problems by running checkGAMSSolvers('MIP')</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >solvers = checkGAMSSolvers(</span><span style="color: rgb(170, 4, 249);">'MIP'</span><span >); disp(solvers);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="EF50CD33" data-testid="output_3" data-width="428" data-height="157" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">  Columns 1 through 12

    'AMPL'    'BARON'    'BDMLP'    'BENCH'    'CBC'    'CONVERT'    'CPLEX'    'EXAMINER'    'GAMSCHK'    'GUROBI'    'KESTREL'    'LINDO'

  Columns 13 through 22

    'LINDOGLOBAL'    'LINGO'    'LOCALSOLVER'    'MOSEK'    'MPECDUMP'    'OsiCplex'    'OsiGurobi'    'OsiMosek'    'OsiXpress'    'PYOMO'

  Columns 23 through 25

    'SCIP'    'XA'    'XPRESS'</div></div></div></div></div><div  class = 'S10'><span>We then run findMustLWithGAMS.m and findMustUWithGAMS.m. </span></div><div  class = 'S2'><span style=' font-weight: bold;'>i) MustL Set: </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[mustLSet, pos_mustL] = findMustLWithGAMS(model, minFluxes_WT, maxFluxes_WT, </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsFindMustL'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'MustL'</span><span >, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div></div><div  class = 'S10'><span>Note that the folder "TestOptForce" was created. Inside this folder, two additional folders were created: "InputsMustL" and "OutputsFindMustL". In the inputs folder you will find all the inputs required to run the GAMS function "findMustL.gms". Additionally, in the outputs folder you will find the mustL set found, which were saved in two files (.xls and .txt), as well as other files generated automatically by GAMS. Furthermore, a report which summarize all the inputs and outputs used during your running was generated. The name of the report will be in this format "report-Day-Month-Year-Hour-Minutes". So, you can mantain a chronological order of your experiments</span></div><div  class = 'S2'><span>We display the reactions that belongs to the mustL set</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >disp(mustLSet)</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B299BE56" data-testid="output_4" data-width="428" data-height="409" data-hashorizontaloverflow="false" style="width: 458px; max-height: 420px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'R11'
    'R26'
    'R37'
    'R38'
    'R39'
    'R40'
    'R41'
    'R42'
    'R43'
    'R46'
    'R48'
    'R49'
    'R50'
    'R51'
    'R52'
    'R53'
    'R54'
    'R55'
    'R56'
    'R57'
    'R58'
    'R59'
    'R60'
    'R61'
    'R73'
    'R74'
    'PSEUDOpyr_1'
    'PSEUDOpep_1'
    'PSEUDOco2_1'</div></div></div></div><div class="inlineWrapper"><div  class = 'S16'></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>ii) MustU set: </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[mustUSet, pos_mustU] = findMustUWithGAMS(model, minFluxes_WT, maxFluxes_WT, </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsFindMustU'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'MustU'</span><span >, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div></div><div  class = 'S10'><span>Note that the folders "InputsMustU" and "OutputsFindMustU" were created. These folders contain the inputs and outputs of findMustUWithGAMS.m, respectively. </span></div><div  class = 'S2'><span>We display the reactions that belongs to the mustU set</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >disp(mustUSet)</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A8EA403E" data-testid="output_5" data-width="428" data-height="185" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'R21'
    'R22'
    'R23'
    'R24'
    'R33'
    'R34'
    'R35'
    'R36'
    'R69'
    'EX_pdo'
    'EX_nh3'
    'EX_so4'
    'SUCt'</div></div></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>B) Finding second order must sets </span></div><div  class = 'S2'><span>First, we define the reactions that will be excluded from the analysis. It it suggested to eliminate reactions found in the previous step as well as exchange reactions</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >constrOpt = struct(</span><span style="color: rgb(170, 4, 249);">'rxnList'</span><span >, {{</span><span style="color: rgb(170, 4, 249);">'EX_gluc'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R75'</span><span >,</span><span style="color: rgb(170, 4, 249);">'EX_suc'</span><span >}}, </span><span style="color: rgb(170, 4, 249);">'values'</span><span >, [-100, 0, 155.5]');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >exchangeRxns = model.rxns(cellfun(@isempty, strfind(model.rxns, </span><span style="color: rgb(170, 4, 249);">'EX_'</span><span >)) == 0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >excludedRxns = unique([mustUSet; mustLSet; exchangeRxns]);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >mustSetFirstOrder = unique([mustUSet; mustLSet]);</span></span></div></div></div><div  class = 'S10'><span>Now, we run the functions for finding second order must sets</span></div><div  class = 'S2'><span style=' font-weight: bold;'>i) MustUU: </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[mustUU, pos_mustUU, mustUU_linear, pos_mustUU_linear] = findMustUUWithGAMS(model, minFluxes_WT, maxFluxes_WT, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >     </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(170, 4, 249);">'excludedRxns'</span><span >, excludedRxns, </span><span style="color: rgb(170, 4, 249);">'mustSetFirstOrder'</span><span >, mustSetFirstOrder, </span><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsFindMustUU'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'MustUU' </span><span >, </span><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div></div><div  class = 'S10'><span>Note that the folders "InputsMustUU" and "OutputsFindMustUU" were created. These folders contain the inputs and outputs of findMustUUWithGAMS.m, respectively. </span></div><div  class = 'S2'><span>We display the reactions that belongs to the mustUU set</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >disp(mustUU);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="938A2A07" data-testid="output_6" data-width="428" data-height="31" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'R30'    'R65'
    'R65'    'R31'</div></div></div></div><div class="inlineWrapper"><div  class = 'S16'></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>ii) MustLL: </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[mustLL, pos_mustLL, mustLL_linear, pos_mustLL_linear] = findMustLLWithGAMS(model, minFluxes_WT, maxFluxes_WT, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >     </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(170, 4, 249);">'excludedRxns'</span><span >, excludedRxns, </span><span style="color: rgb(170, 4, 249);">'mustSetFirstOrder'</span><span >, mustSetFirstOrder, </span><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsFindMustLL'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'MustLL' </span><span >, </span><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div></div><div  class = 'S10'><span>Note that the folders "InputsMustLL" and "OutputsFindMustLL" were created. These folders contain the inputs and outputs of findMustLLWithGAMS.m, respectively. </span></div><div  class = 'S2'><span>We display the reactions that belongs to the mustLL set. In this case, MustLL is an empty array because no reaction was found in the mustLL set.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: normal"><span >disp(mustLL);</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>iii) MustUL: </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[mustUL, pos_mustUL, mustUL_linear, pos_mustUL_linear] = findMustULWithGAMS(model, minFluxes_WT, maxFluxes_WT, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >     </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(170, 4, 249);">'excludedRxns'</span><span >, excludedRxns, </span><span style="color: rgb(170, 4, 249);">'mustSetFirstOrder'</span><span >, mustSetFirstOrder, </span><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsFindMustUL'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'MustUL' </span><span >, </span><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div></div><div  class = 'S10'><span>Note that the folders "InputsMustUL" and "OutputsFindMustUL" were created. These folders contain the inputs and outputs of findMustULWithGAMS.m, respectively. </span></div><div  class = 'S2'><span>We display the reactions that belongs to the mustUL set. In this case, MustUL is an empty array because no reaction was found in the mustUL set.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: normal"><span >disp(mustUL);</span></span></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>TROUBLESHOOTING 1: </span><span> "I didn't find any reaction in my must sets"</span></div><div  class = 'S2'><span style=' font-weight: bold;'>TROUBLESHOOTING 2: </span><span> "I got an error when running the findMustXWithGAMS.m functions (X = L or U or LL or UL or UU depending on the case)"</span></div><h2  class = 'S3'><span>Step 5: OptForce</span></h2><div  class = 'S2'><span style=' font-weight: bold;'>TIMING: </span><span>This task should take from a few seconds to a few hours depending on the size of your reconstruction</span></div><div  class = 'S2'><span>We define constraints and we define "K" the number of interventions allowed, "nSets" the maximum number of sets to find, and "targetRxn" the reaction producing the metabolite of interest (in this case, succinate). </span></div><div  class = 'S2'><span>Additionally, we define the mustU set as the union of the reactions that must be upregulated in both first and second order must sets; and mustL set as the union of the reactions that must be downregulated in both first and second order must sets </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >mustU = unique(union(mustUSet, mustUU));</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >mustL = unique(union(mustLSet, mustLL));</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >targetRxn = </span><span style="color: rgb(170, 4, 249);">'EX_suc'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >k = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >nSets = 1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >constrOpt = struct(</span><span style="color: rgb(170, 4, 249);">'rxnList'</span><span >, {{</span><span style="color: rgb(170, 4, 249);">'EX_gluc'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R75'</span><span >}}, </span><span style="color: rgb(170, 4, 249);">'values'</span><span >, [-100, 0]);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >[optForceSets, posOptForceSets, typeRegOptForceSets] = optForceWithGAMS(model, targetRxn, mustU, mustL, minFluxes_WT, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    maxFluxes_WT, minFluxes_MT, maxFluxes_MT, </span><span style="color: rgb(170, 4, 249);">'k'</span><span >, k, </span><span style="color: rgb(170, 4, 249);">'nSets'</span><span >, nSets, </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsOptForce'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OptForce'</span><span >, </span><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1,</span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div></div><div  class = 'S10'><span>Note that the folders "InputsOptForce" and "OutputsOptForce" were created. These folders contain the inputs and outputs of optForceWithGAMS.m, respectively.</span></div><div  class = 'S2'><span>We display the reactions found by optForce</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >disp(optForceSets)</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="1AA94EDD" data-testid="output_7" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'SUCt'</div></div></div></div><div class="inlineWrapper"><div  class = 'S16'></div></div></div><div  class = 'S10'><span>The reaction found was "SUCt", i.e. a transporter for succinate (a very intuitive solution).</span></div><div  class = 'S2'><span>Next, we will increase "k" and we will exclude "SUCt" from upregulations to found non-intuitive solutions. </span></div><div  class = 'S2'><span style=' font-weight: bold;'>TIP: </span><span>Sometimes the product is at the end of a long linear pathway. In that case, the recomendation is to also exclude most reactions on the linear pathway. Essential reactions and reactions not associated with any gene should also be excluded. </span></div><div  class = 'S2'><span>We will only search for the 20 best solutions, but you can try with a higher number.</span></div><div  class = 'S2'><span>We will change the runID to save both resutls (k = 1 and K = 2) in diffetent folders</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >k = 2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >nSets = 20;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >runID = </span><span style="color: rgb(170, 4, 249);">'TestOptForce2'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >excludedRxns = struct(</span><span style="color: rgb(170, 4, 249);">'rxnList'</span><span >, {{</span><span style="color: rgb(170, 4, 249);">'SUCt'</span><span >}}, </span><span style="color: rgb(170, 4, 249);">'typeReg'</span><span >,</span><span style="color: rgb(170, 4, 249);">'U'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >[optForceSets, posOptForceSets, typeRegOptForceSets] = optForceWithGAMS(model, targetRxn, mustU, mustL, minFluxes_WT, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    maxFluxes_WT, minFluxes_MT, maxFluxes_MT, </span><span style="color: rgb(170, 4, 249);">'k'</span><span >, k, </span><span style="color: rgb(170, 4, 249);">'nSets'</span><span >, nSets, </span><span style="color: rgb(170, 4, 249);">'constrOpt'</span><span >, constrOpt, </span><span style="color: rgb(170, 4, 249);">'excludedRxns'</span><span >, excludedRxns, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'runID'</span><span >, runID, </span><span style="color: rgb(170, 4, 249);">'outputFolder'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OutputsOptForce'</span><span >, </span><span style="color: rgb(170, 4, 249);">'outputFileName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'OptForce'</span><span >, </span><span style="color: rgb(170, 4, 249);">'solverName'</span><span >, </span><span style="color: rgb(170, 4, 249);">'cplex'</span><span >, </span><span style="color: rgb(170, 4, 249);">'printExcel'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'printText'</span><span >, 1, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'printReport'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepInputs'</span><span >, 1, </span><span style="color: rgb(170, 4, 249);">'keepGamsOutputs'</span><span >, 1,</span><span style="color: rgb(170, 4, 249);">'verbose'</span><span >, 0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div></div><div  class = 'S10'><span>Note that the folders "InputsOptForce" and "OutputsOptForce" were created inside TestOptForce2. These folders contain the inputs and outputs of optForceWithGAMS.m, respectively.</span></div><div  class = 'S2'><span>We display the reactions found by optForce</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >disp(optForceSets)</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A836E66B" data-testid="output_8" data-width="428" data-height="283" data-hashorizontaloverflow="false" style="width: 458px; max-height: 294px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'R23'    'R26'
    'R24'    'R25'
    'R24'    'R26'
    'R23'    'R25'
    'R23'    'R26'
    'R23'    'R63'
    'R22'    'R25'
    'R24'    'R63'
    'R21'    'R26'
    'R24'    'R26'
    'R22'    'R26'
    'R22'    'R63'
    'R22'    'R26'
    'R21'    'R63'
    'R21'    'R25'
    'R21'    'R26'
    'R24'    'R4' 
    'R23'    'R4' 
    'R22'    'R4' 
    'R21'    'R4' </div></div></div></div></div><h2  class = 'S3'><span>TIMING</span></h2><ol  class = 'S4'><li  class = 'S5'><span>STEP 1: &lt; 1 second</span></li><li  class = 'S5'><span>STEP 2: ~ 1-2 seconds</span></li><li  class = 'S5'><span>STEP 3:  ~ 2-5 seconds</span></li><li  class = 'S5'><span>STEP 4: ~ 10-20 seconds</span></li><li  class = 'S5'><span>STEP 5: ~ 10-20 seconds</span></li></ol><h2  class = 'S3'><span>TROUBLESHOOTING</span></h2><div  class = 'S2'><span>1) problem: "I didn't find any reaction in my must sets"</span></div><div  class = 'S2'><span>possible reason: the wild-type or mutant strain is not constrained enough. </span></div><div  class = 'S2'><span>solution: add more constraints to your strains until you find differences in your reaction ranges. If you don't find any differences, it is better to change the approach and use another algorithm. Also, it is possible that you won't find second order must set (like in this tutorial). You can check if the algorithm is working by defining excludedRxns as an empty array. If the algorihm is working well, you should see how the must sets are found in each iteration.</span></div><div  class = 'S2'><span></span></div><div  class = 'S2'><span>2) problem: "I got an error when running the findMust functions"</span></div><div  class = 'S2'><span>possible reason: inputs are not defined well or GAMS was not installed correctly</span></div><div  class = 'S2'><span>solution: verify your inputs. To verify that GAMS was installed correctly, in MATLAB, use the function which('gams'). If a folder is shown, GAMS was installed correctly.</span></div><h2  class = 'S3'><span>ANTICIPATED RESULTS</span></h2><div  class = 'S2'><span>In this tutorial some folders will be created inside the folder called "runID" to store inputs and outputs of the optForce functions (findMustUWithGams.m, findMustLWithGams.m, findMustUUWithGams.m, findMustLLWithGams.m, findMustULWithGams.m, optForceWithGams.m)</span></div><div  class = 'S2'><span>In this case runID = 'TestOptForce', so inside this folder the following folders will be created:</span></div><div  class = 'S2'><span style=' font-family: monospace;'>├── CurrentFolder</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   ├── TestOptForce</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── InputsFindMustL</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── OutputsFindMustL</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── InputsFindMustU</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── OutputsFindMustU</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── InputsFindMustLL</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── OutputsFindMustLL</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── InputsFindMustUU</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── OutputsFindMustUU</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── InputsFindMustUL</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── OutputsFindMustUL</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   ├── InputsOptForce</span></div><div  class = 'S2'><span style=' font-family: monospace;'>|   |   └── OutputsOptForce</span></div><div  class = 'S2'><span style=' font-family: monospace;'></span></div><div  class = 'S2'><span>The input folders contain inputs (.txt and .gdx files) for running the GAMS functions to solve each one of the bilevel problems. Output folders contain results of the algorithms (.xls and .txt files) as well as a report (.txt) summarizing the outcomes of the steps performed during the execution of the optForce functions. Additionally, some auto-generated files will be stored (.gdx and .lst files). These files are generated by GAMS automatically and contain values for variables and the summary of the GAMS execution, respectively.</span></div><div  class = 'S2'><span>The optForce algorithm will find sets of reactions that should increase the production of your target. The first sets found should be the best ones because the production rate will be the highest. The last ones should be the worse because the production rete will be slower. Be aware that some sets could not guarante a minimum production rate for your target, so you always have to check the minimum production rate. You can do this using the function testOptForceSol.m. Some sets could allow a higher growth rate than others, so keep in mind this too when deciding which set is better.</span></div><h2  class = 'S3'><span>Acknowledgments</span></h2><div  class = 'S2'><span>I would to thanks to the research group of Costas D. Maranas who provided the GAMS functions to solve this example. In particular I would like to thank to Chiam Yu Ng who kindly provided examples for using GAMS and support for writing GAMS functions.</span></div><h2  class = 'S3'><span>References</span></h2><div  class = 'S2'><span>[1] Ranganathan S, Suthers PF, Maranas CD (2010) OptForce: An Optimization Procedure for Identifying All Genetic Manipulations Leading to Targeted Overproductions. PLOS Computational Biology 6(4): e1000744. https://doi.org/10.1371/journal.pcbi.1000744.</span></div><div  class = 'S2'><span>[2] Maciek R. Antoniewicz, David F. Kraynie, Lisa A. Laffend, Joanna González-Lergier, Joanne K. Kelleher, Gregory Stephanopoulos, Metabolic flux analysis in a nonstationary system: Fed-batch fermentation of a high yielding strain of E. coli producing 1,3-propanediol, Metabolic Engineering, Volume 9, Issue 3, May 2007, Pages 277-292, ISSN 1096-7176, https://doi.org/10.1016/j.ymben.2007.01.003.</span></div><div  class = 'S2'><span></span></div><div  class = 'S2'></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% OptForce GAMS
%% Author: Sebastián N. Mendoza,  Center for Mathematical Modeling, University of Chile. snmendoz@uc.cl
%% *Reviewers(s): Chiam Yu Ng (Costas D. Maranas group),* Lin Wang *(Costas D. Maranas group)*
%% *INTRODUCTION:*
% In this tutorial we will run optForce. For a detailed description of the procedure, 
% please see [1]. Briefly, the problem is to find a set of interventions of size 
% "K" such that when these interventions are applied to a wild-type strain, the 
% mutant created will produce a particular target of interest in a higher rate 
% than the wild-type strain. The interventions could be knockouts (lead to zero 
% the flux for a particular reaction), upregulations (increase the flux for a 
% particular reaction) and downregulations (decrease the flux for a particular 
% reaction). 
% 
% For example, imagine that we would like to increase the production of succinate 
% in Escherichia coli. Which are the interventions needed to increase the production 
% of succinate? We will approach this problem in this tutorial and we will see 
% how each of the steps of OptForce are solved. 
%% MATERIALS
%% EQUIPMENT
%% 
% # MATLAB
% # GAMS
% # A solver for Mixed Integer Linear Programming (MILP) problems. For example, 
% Gurobi.
%% *EQUIPMENT SETUP*
% *GAMS*: Install the latest version of GAMS available. This is important as 
% some GAMS functions are not supported in versions older than  24.7.1. During 
% GAMS installation make sure to select the option "Add GAMS directory to PATH 
% environment variable". In MATLAB, use functions addpath() and savepath() for 
% adding and saving the directory where GAMS was installed. For a more detailed 
% description of GAMS installation and setup, you can watch the official tutorial 
% of GDXMRW at https://www.youtube.com/watch?v=HPn_q8nlktE
%% PROCEDURE
% The proceduce consists on the following steps
% 
% 1) Maximize specific growth rate and product formation.
% 
% 2) Define constraints for both wild-type and mutant strain: 
% 
% 3) Perform flux variability analysis for both wild-type and mutant strain.
% 
% 4) Find must sets,  i.e, reactions that MUST increase or decrease their flux 
% in order to achieve the phenotype in the mutant strain. 
%% Figure 1.
% 
% 
% 5) Find the interventions needed that will ensure a increased production of 
% the target of interest
% 
% Now, we will approach each step in detail.
%% STEP 1: Maximize specific growth rate and product formation
% First, we load the model. This model comprises only 90 reactions, which describe 
% the central metabolism of E. coli [2].
% 
% Then, we change the objective function to maximize biomass ("R75"). We also 
% change the lower bounds, so E. coli will be able to consume glucose, oxygen, 
% sulfate, ammomium, citrate and glycerol.

global TUTORIAL_INIT_CB;
if ~isempty(TUTORIAL_INIT_CB) && TUTORIAL_INIT_CB==1
    initCobraToolbox(false) % false, as we don't want to update
    changeCobraSolver('gurobi','all');
end

pathTutorial = which('tutorial_OptForceGAMS.mlx');
pathstr = fileparts(pathTutorial);
cd(pathstr)

modelFileName = 'AntCore.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
model = readCbModel(modelFileName);

model.c(strcmp(model.rxns,'R75')) = 1;
model = changeRxnBounds(model, 'EX_gluc', -100, 'l'); 
model = changeRxnBounds(model, 'EX_o2', -100, 'l'); 
model = changeRxnBounds(model, 'EX_so4', -100, 'l'); 
model = changeRxnBounds(model, 'EX_nh3', -100, 'l'); 
model = changeRxnBounds(model, 'EX_cit', -100, 'l'); 
model = changeRxnBounds(model, 'EX_glyc', -100, 'l'); 
%% 
% Then, we calculate the maximum specific growth rate and the maximum production 
% rate for succinate

growthRate = optimizeCbModel(model); 
fprintf('The maximum growth rate is %1.2f', growthRate.f);

model = changeObjective(model, 'EX_suc');
maxSucc = optimizeCbModel(model);
fprintf('The maximum production rate of succinate is %1.2f', maxSucc.f);
%% 
% *TIP:* The biomass reaction is usually set to 1%-10% of maximum theoretical 
% biomass yield when running the following steps, to prevent solutions with not 
% biomass formation
%% 
% # maximizing product formation
% # finding MUST sets of second order
% # finding FORCE sets
%% STEP 2: Define constraints for both wild-type and mutant strain
% *TIMING*: This step should take a few days or weeks, depending on the information 
% available for your species. 
% 
% *CRITICAL STEP*: This is a manual task, so you should search for information 
% in articles or even perform your own experiments. You can also make assumptions 
% for describing the phenotypes of both strains which will make this task a little 
% faster but make sure to have two strains different enough, because you should 
% be able to find differences in reactions ranges. 
% 
% First, we load the model. This model comprises only 90 reactions, which describe 
% the central metabolism of E. coli [2].
% 
% Then, we change the objective function to maximize biomass ("R75"). We also 
% change the lower bounds, so E. coli will be able to consume glucose, oxygen, 
% sulfate, ammomium, citrate and glycerol. 
% 
% We define constraints for each strain

Constr_WT = struct('rxnList', {{'R75'}}, 'rxnValues', 14, 'rxnBoundType', 'b');
Constr_MT = struct('rxnList', {{'R75','EX_suc'}}, 'rxnValues', [0, 155.55], 'rxnBoundType', 'bb');
%% Step 3: Flux Variability Analysis
% *TIMING*: This task should take from a few seconds to a few hours depending 
% on the size of your reconstruction
% 
% We  run the FVA analysis for both strains

[minFluxes_WT, maxFluxes_WT, minFluxes_MT, maxFluxes_MT,~,~] = FVAOptForce(model, Constr_WT, Constr_MT);
disp([minFluxes_WT, maxFluxes_WT, minFluxes_MT, maxFluxes_MT]);
%% 
% Now, the run the second step of OptForce.
%% Step 4: Find Must Sets
% *TIMING:* This task should take from a few seconds to a few hours depending 
% on the size of your reconstruction
% 
% First, we define an ID for this run. Each time you run the functions associated 
% to the optForce procedure, some folders can be generated to store inputs used 
% in that run. Outputs are stored as well. These folder will be located inside 
% the folder defined by your run ID. Thus, if your runID is ''TestOptForce", the 
% structure of the folders will be the following:
% 
% |├── CurrentFolder|
% 
% ||   ├── TestOptForce|
% 
% ||   |   ├── Inputs|
% 
% ||   |   └── Outputs|
% 
% To avoid the generation of inputs and outputs folders, set keepInputs = 0, 
% printExcel = 0, printText = 0 and keepGamsOutputs = 0
% 
% Also, a report of the run is generated each time you run the functions associated 
% to the optForce procedure. So, the idea is to give a different runID each time 
% you run the functions, so you will be able to see the report (inputs used, outputs 
% generated, errors in the run) for each run.
% 
% We define then our runID

runID = 'TestOptForce';
%% 
% Fow now, only functions to find first and second order must sets are supported. 
% As depicted in *Figure 1*, the first order must sets are MUSTU and MUSTL; and 
% second order must sets are MUSTUU, MUSTLL and MUSTUL
% 
% *A) Finding first order must sets*
% 
% We define constraints

constrOpt = struct('rxnList', {{'EX_gluc','R75','EX_suc'}}, 'values', [-100, 0, 155.5]');
%% 
% We then run the functions findMustLWithGAMS.m and findMustUWithGAMS.m that 
% will find mustU and mustL sets, respectively.
% 
% Important: To run these function you will need a solver able to solve Mixed 
% Integer Linear Programming (MILP or MIP) problems. Some popular options are: 
% cplex and gurobi. You can see which gams solvers are available in your systems 
% to solve MIP problems by running checkGAMSSolvers('MIP')

solvers = checkGAMSSolvers('MIP'); disp(solvers);
%% 
% We then run findMustLWithGAMS.m and findMustUWithGAMS.m. 
% 
% *i) MustL Set:* 

[mustLSet, pos_mustL] = findMustLWithGAMS(model, minFluxes_WT, maxFluxes_WT, 'constrOpt', constrOpt, ...
'solverName', 'cplex', 'runID', runID, 'outputFolder', 'OutputsFindMustL', 'outputFileName', 'MustL', ...
'printExcel', 1, 'printText', 1, 'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1, 'verbose', 0);
%% 
% Note that the folder "TestOptForce" was created. Inside this folder, two additional 
% folders were created: "InputsMustL" and "OutputsFindMustL". In the inputs folder 
% you will find all the inputs required to run the GAMS function "findMustL.gms". 
% Additionally, in the outputs folder you will find the mustL set found, which 
% were saved in two files (.xls and .txt), as well as other files generated automatically 
% by GAMS. Furthermore, a report which summarize all the inputs and outputs used 
% during your running was generated. The name of the report will be in this format 
% "report-Day-Month-Year-Hour-Minutes". So, you can mantain a chronological order 
% of your experiments
% 
% We display the reactions that belongs to the mustL set

disp(mustLSet)

%% 
% *ii) MustU set:* 

[mustUSet, pos_mustU] = findMustUWithGAMS(model, minFluxes_WT, maxFluxes_WT, 'constrOpt', constrOpt, ...
'solverName', 'cplex', 'runID', runID, 'outputFolder', 'OutputsFindMustU', 'outputFileName', 'MustU', ...
'printExcel', 1, 'printText', 1, 'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1, 'verbose', 0);
%% 
% Note that the folders "InputsMustU" and "OutputsFindMustU" were created. These 
% folders contain the inputs and outputs of findMustUWithGAMS.m, respectively. 
% 
% We display the reactions that belongs to the mustU set

disp(mustUSet)
%% 
% *B) Finding second order must sets* 
% 
% First, we define the reactions that will be excluded from the analysis. It 
% it suggested to eliminate reactions found in the previous step as well as exchange 
% reactions

constrOpt = struct('rxnList', {{'EX_gluc','R75','EX_suc'}}, 'values', [-100, 0, 155.5]');
exchangeRxns = model.rxns(cellfun(@isempty, strfind(model.rxns, 'EX_')) == 0);
excludedRxns = unique([mustUSet; mustLSet; exchangeRxns]);
mustSetFirstOrder = unique([mustUSet; mustLSet]);
%% 
% Now, we run the functions for finding second order must sets
% 
% *i) MustUU:* 

[mustUU, pos_mustUU, mustUU_linear, pos_mustUU_linear] = findMustUUWithGAMS(model, minFluxes_WT, maxFluxes_WT, ...
     'constrOpt', constrOpt, 'excludedRxns', excludedRxns, 'mustSetFirstOrder', mustSetFirstOrder, 'solverName', 'cplex',...
    'runID', runID, 'outputFolder', 'OutputsFindMustUU', 'outputFileName', 'MustUU' , 'printExcel', 1, 'printText', 1, ...
    'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1, 'verbose', 0);
%% 
% Note that the folders "InputsMustUU" and "OutputsFindMustUU" were created. 
% These folders contain the inputs and outputs of findMustUUWithGAMS.m, respectively. 
% 
% We display the reactions that belongs to the mustUU set

disp(mustUU);

%% 
% *ii) MustLL:* 

[mustLL, pos_mustLL, mustLL_linear, pos_mustLL_linear] = findMustLLWithGAMS(model, minFluxes_WT, maxFluxes_WT, ...
     'constrOpt', constrOpt, 'excludedRxns', excludedRxns, 'mustSetFirstOrder', mustSetFirstOrder, 'solverName', 'cplex',...
    'runID', runID, 'outputFolder', 'OutputsFindMustLL', 'outputFileName', 'MustLL' , 'printExcel', 1, 'printText', 1, ...
    'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1, 'verbose', 0);
%% 
% Note that the folders "InputsMustLL" and "OutputsFindMustLL" were created. 
% These folders contain the inputs and outputs of findMustLLWithGAMS.m, respectively. 
% 
% We display the reactions that belongs to the mustLL set. In this case, MustLL 
% is an empty array because no reaction was found in the mustLL set.

disp(mustLL);
%% 
% *iii) MustUL:* 

[mustUL, pos_mustUL, mustUL_linear, pos_mustUL_linear] = findMustULWithGAMS(model, minFluxes_WT, maxFluxes_WT, ...
     'constrOpt', constrOpt, 'excludedRxns', excludedRxns, 'mustSetFirstOrder', mustSetFirstOrder, 'solverName', 'cplex',...
    'runID', runID, 'outputFolder', 'OutputsFindMustUL', 'outputFileName', 'MustUL' , 'printExcel', 1, 'printText', 1, ...
    'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1, 'verbose', 0);
%% 
% Note that the folders "InputsMustUL" and "OutputsFindMustUL" were created. 
% These folders contain the inputs and outputs of findMustULWithGAMS.m, respectively. 
% 
% We display the reactions that belongs to the mustUL set. In this case, MustUL 
% is an empty array because no reaction was found in the mustUL set.

disp(mustUL);
%% 
% *TROUBLESHOOTING 1:*  "I didn't find any reaction in my must sets"
% 
% *TROUBLESHOOTING 2:*  "I got an error when running the findMustXWithGAMS.m 
% functions (X = L or U or LL or UL or UU depending on the case)"
%% Step 5: OptForce
% *TIMING:* This task should take from a few seconds to a few hours depending 
% on the size of your reconstruction
% 
% We define constraints and we define "K" the number of interventions allowed, 
% "nSets" the maximum number of sets to find, and "targetRxn" the reaction producing 
% the metabolite of interest (in this case, succinate). 
% 
% Additionally, we define the mustU set as the union of the reactions that must 
% be upregulated in both first and second order must sets; and mustL set as the 
% union of the reactions that must be downregulated in both first and second order 
% must sets 

mustU = unique(union(mustUSet, mustUU));
mustL = unique(union(mustLSet, mustLL));
targetRxn = 'EX_suc';
k = 1;
nSets = 1;
constrOpt = struct('rxnList', {{'EX_gluc','R75'}}, 'values', [-100, 0]);

[optForceSets, posOptForceSets, typeRegOptForceSets] = optForceWithGAMS(model, targetRxn, mustU, mustL, minFluxes_WT, ...
    maxFluxes_WT, minFluxes_MT, maxFluxes_MT, 'k', k, 'nSets', nSets, 'constrOpt', constrOpt, ...
    'runID', runID, 'outputFolder', 'OutputsOptForce', 'outputFileName', 'OptForce', 'solverName', 'cplex', 'printExcel', 1, 'printText', 1, ...
    'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1,'verbose', 0);

%% 
% Note that the folders "InputsOptForce" and "OutputsOptForce" were created. 
% These folders contain the inputs and outputs of optForceWithGAMS.m, respectively.
% 
% We display the reactions found by optForce

disp(optForceSets)

%% 
% The reaction found was "SUCt", i.e. a transporter for succinate (a very intuitive 
% solution).
% 
% Next, we will increase "k" and we will exclude "SUCt" from upregulations to 
% found non-intuitive solutions. 
% 
% *TIP:* Sometimes the product is at the end of a long linear pathway. In that 
% case, the recomendation is to also exclude most reactions on the linear pathway. 
% Essential reactions and reactions not associated with any gene should also be 
% excluded. 
% 
% We will only search for the 20 best solutions, but you can try with a higher 
% number.
% 
% We will change the runID to save both resutls (k = 1 and K = 2) in diffetent 
% folders

k = 2;
nSets = 20;
runID = 'TestOptForce2';
excludedRxns = struct('rxnList', {{'SUCt'}}, 'typeReg','U');
[optForceSets, posOptForceSets, typeRegOptForceSets] = optForceWithGAMS(model, targetRxn, mustU, mustL, minFluxes_WT, ...
    maxFluxes_WT, minFluxes_MT, maxFluxes_MT, 'k', k, 'nSets', nSets, 'constrOpt', constrOpt, 'excludedRxns', excludedRxns, ...
    'runID', runID, 'outputFolder', 'OutputsOptForce', 'outputFileName', 'OptForce', 'solverName', 'cplex', 'printExcel', 1, 'printText', 1, ...
    'printReport', 1, 'keepInputs', 1, 'keepGamsOutputs', 1,'verbose', 0);

%% 
% Note that the folders "InputsOptForce" and "OutputsOptForce" were created 
% inside TestOptForce2. These folders contain the inputs and outputs of optForceWithGAMS.m, 
% respectively.
% 
% We display the reactions found by optForce

disp(optForceSets)
%% TIMING
%% 
% # STEP 1: < 1 second
% # STEP 2: ~ 1-2 seconds
% # STEP 3:  ~ 2-5 seconds
% # STEP 4: ~ 10-20 seconds
% # STEP 5: ~ 10-20 seconds
%% TROUBLESHOOTING
% 1) problem: "I didn't find any reaction in my must sets"
% 
% possible reason: the wild-type or mutant strain is not constrained enough. 
% 
% solution: add more constraints to your strains until you find differences 
% in your reaction ranges. If you don't find any differences, it is better to 
% change the approach and use another algorithm. Also, it is possible that you 
% won't find second order must set (like in this tutorial). You can check if the 
% algorithm is working by defining excludedRxns as an empty array. If the algorihm 
% is working well, you should see how the must sets are found in each iteration.
% 
% 
% 
% 2) problem: "I got an error when running the findMust functions"
% 
% possible reason: inputs are not defined well or GAMS was not installed correctly
% 
% solution: verify your inputs. To verify that GAMS was installed correctly, 
% in MATLAB, use the function which('gams'). If a folder is shown, GAMS was installed 
% correctly.
%% ANTICIPATED RESULTS
% In this tutorial some folders will be created inside the folder called "runID" 
% to store inputs and outputs of the optForce functions (findMustUWithGams.m, 
% findMustLWithGams.m, findMustUUWithGams.m, findMustLLWithGams.m, findMustULWithGams.m, 
% optForceWithGams.m)
% 
% In this case runID = 'TestOptForce', so inside this folder the following folders 
% will be created:
% 
% |├── CurrentFolder|
% 
% ||   ├── TestOptForce|
% 
% ||   |   ├── InputsFindMustL|
% 
% ||   |   ├── OutputsFindMustL|
% 
% ||   |   ├── InputsFindMustU|
% 
% ||   |   ├── OutputsFindMustU|
% 
% ||   |   ├── InputsFindMustLL|
% 
% ||   |   ├── OutputsFindMustLL|
% 
% ||   |   ├── InputsFindMustUU|
% 
% ||   |   ├── OutputsFindMustUU|
% 
% ||   |   ├── InputsFindMustUL|
% 
% ||   |   ├── OutputsFindMustUL|
% 
% ||   |   ├── InputsOptForce|
% 
% ||   |   └── OutputsOptForce|
% 
% 
% 
% The input folders contain inputs (.txt and .gdx files) for running the GAMS 
% functions to solve each one of the bilevel problems. Output folders contain 
% results of the algorithms (.xls and .txt files) as well as a report (.txt) summarizing 
% the outcomes of the steps performed during the execution of the optForce functions. 
% Additionally, some auto-generated files will be stored (.gdx and .lst files). 
% These files are generated by GAMS automatically and contain values for variables 
% and the summary of the GAMS execution, respectively.
% 
% The optForce algorithm will find sets of reactions that should increase the 
% production of your target. The first sets found should be the best ones because 
% the production rate will be the highest. The last ones should be the worse because 
% the production rete will be slower. Be aware that some sets could not guarante 
% a minimum production rate for your target, so you always have to check the minimum 
% production rate. You can do this using the function testOptForceSol.m. Some 
% sets could allow a higher growth rate than others, so keep in mind this too 
% when deciding which set is better.
%% Acknowledgments
% I would to thanks to the research group of Costas D. Maranas who provided 
% the GAMS functions to solve this example. In particular I would like to thank 
% to Chiam Yu Ng who kindly provided examples for using GAMS and support for writing 
% GAMS functions.
%% References
% [1] Ranganathan S, Suthers PF, Maranas CD (2010) OptForce: An Optimization 
% Procedure for Identifying All Genetic Manipulations Leading to Targeted Overproductions. 
% PLOS Computational Biology 6(4): e1000744. https://doi.org/10.1371/journal.pcbi.1000744.
% 
% [2] Maciek R. Antoniewicz, David F. Kraynie, Lisa A. Laffend, Joanna González-Lergier, 
% Joanne K. Kelleher, Gregory Stephanopoulos, Metabolic flux analysis in a nonstationary 
% system: Fed-batch fermentation of a high yielding strain of E. coli producing 
% 1,3-propanediol, Metabolic Engineering, Volume 9, Issue 3, May 2007, Pages 277-292, 
% ISSN 1096-7176, https://doi.org/10.1016/j.ymben.2007.01.003.
% 
% 
% 
%
##### SOURCE END #####
-->
</div></body></html>